Anatomy of the binary black hole recoil: A multipolar analysis 
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We present a multipolar analysis of the recoil velocity computed in recent numerical simulations 
of binary black hole coalescence, for both unequal masses and non-zero, non-precessing spins. We 
show that multipole moments up to and including l = 4 are sufficient to accurately reproduce the 
final recoil velocity (~ 98%) and that only a few dominant modes contribute significantly to it 
(— 95%). We describe how the relative amplitude, and more importantly, the relative phase, of 
these few modes control the way in which the recoil builds up throughout the inspiral, merger, and 
ring-down phases. We also find that the numerical results can be reproduced, to a high level of 
accuracy, by an effective Newtonian formula for the multipole moments obtained by replacing in 
the Newtonian formula the radial separation with an effective radius computed from the numerical 
data. Beyond the merger, the numerical results are reproduced by a superposition of three Kerr 
quasi-normal modes. Analytic formulae, obtained by expressing the multipole moments in terms 
of the fundamental QNMs of a Kerr BH, are able to explain the onset and amount of “anti-kick” 
for each of the simulations. Lastly, we apply this multipolar analysis to understand the remarkable 
difference between the amplitudes of planar and non-planar kicks for equal-mass spinning black 
holes. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.70.Bw, 04.25.Nx, 04.30.-w 


I. INTRODUCTION 

After the recent breakthrough in numerical relativity 
(NR) [1-3], a number of different groups [4-6] are now 
able to evolve binary black holes (BHs) through merger. 
Recently, a great deal of effort has been directed towards 
the computation of the recoil velocity of the final BH [7- 
15]. The fundamental cause of this recoil is a net lin- 
ear momentum flux in the gravitational radiation, due to 
some asymmetry in the system [16-18], typically unequal 
masses or spins in the case of BH binaries. The recoil has 
astrophysical importance because it can affect the growth 
of supermassive black holes (SMBHs) through hierarchi- 
cal mergers [19, 20]. In those scenarios dark-matter halos 
grow through hierarchical mergers. The SMBHs at the 
centers of such haloes are thought to merge unless they 
had been kicked out of the gravitational potential well 
because the recoil velocity gained in a previous merger is 
larger than the halo’s escape velocity. 

Other astrophysical implications include the displace- 
ment of the SMBH, along with its gaseous accretion disk, 
forming an “off-center” quasar. These quasars might also 
have emission lines highly red- or blue-shifted relative to 
the host galaxy due to the doppler shift of the recoil 
velocity. Additionally, these displaced SMBHs could in 
turn displace a significant amount of stellar mass from 
the galactic nucleus as they sink back to the center via 
dynamical friction, forming a depleted core of missing 
mass on the order of twice the SMBH mass [20-22]. 

Until now, numerical simulations have computed re- 


coil velocities for non-spinning unequal-mass BH binary 
systems [7-9] in the range m 2 /mi = (1 - - • 4), where 
mi and m 2 are the individual BH masses; for spinning, 
non-precessing binary BHs [10-12], and also for precess- 
ing BHs with both equal [13, 14] as well as unequal 
masses [15]. Quite interestingly, there exist initial spin 
configurations for which the recoil velocity can be quite 
large, e.g., > 2500 km/sec [13-15]. However, it is not yet 
clear whether those very large recoil velocities are astro- 
physically likely [24-26]. So far, due to limited computa- 
tional resources, the numerical simulations have explored 
a rather small portion of the total parameter space. 

Analytic calculations, based on the post-Newtonian 
(PN) expansion of Einstein field equations [27] and its 
siblings [28-33], have made predictions for the recoil ve- 
locity [34-38] before the NR breakthrough. Since the 
majority of the linear momentum flux is emitted during 
the merger and ringdown (RD) phases, it is difficult to 
make definitive predictions for the recoil using only ana- 
lytic methods. These methods need to be somehow cali- 
brated to the NR results, so that they can be accurately 
extended during the transition from inspiral to RD. So 
far, in the non-spinning case, the PN model [37] has pro- 
vided results consistent with NR all along the adiabatic 
inspiral; the effective-one-body (EOB) model [28, 29, 31] 
can reproduce the total recoil, i.e., also the contribution 
from the RD phase, but with large uncertainties [38]. 
More recently, Ref. [23] provided the first estimates of the 
distribution of recoil velocities from spinning BH merg- 
ers using the EOB model, calibrated to the NR results. 
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In Ref. [39], perturbative calculations which make use 
of the so-called close-limit approximation [40] have been 
used to predict the recoil for unequal-mass binary BHs 
moving on circular and eccentric orbits. 

In this paper we present a diagnostic of the physics of 
the recoil, trying to understand how it accumulates dur- 
ing the inspiral, merger, and RD phases. We use several 
numerical simulations of non-spinning, unequal mass bi- 
nary systems, as well as spinning, non-precessing binary 
systems obtained by the Goddard numerical relativity 
group. 

We frame our understanding using the multipolar for- 
malism originally laid out by Thorne [41-45]. We work 
out which (pairs of) multipole moments contribute most 
significantly to the recoil. We employ analytic, but lead- 
ing order, formulae for the multipole moments of the lin- 
ear momentum flux during the inspiral phase, and ex- 
press the multipole moments in terms of a linear su- 
perposition of quasi-normal modes (QNMs) during the 
RD phase [46]. These analysis tools help us understand 
why for some binary mass and spin configurations the 
so-called “anti-kick” is present, and what determines the 
difference between the peak and the final recoil velocity. 
What we learn in this study will be used in a forthcoming 
paper to improve the analytic models, so that they can 
be used to interpolate between NR results, efficiently and 
accurately covering the entire parameter space. 

This paper is organized as follow. In Sec. II, after in- 
troducing our definitions and notations, we review the 
binary parameters used in the numerical simulations and 
examine the main features of the numerical runs. In 
Sec. Ill we discuss the multipolar expansion of the linear 
momentum, angular momentum and energy fluxes given 
in terms of the symmetric trace-free radiative mass and 
current moments, and show how to compute those fluxes 
from the multipole decomposition of the Weyl scalar 
^ 4 . In Sec. IV, we analyse the multipole content of the 
numerical waveforms during the inspiral and ring-down 
phases. In Sec. V we show that, by properly normaliz- 
ing the binary radial separation, the multipole moments 
computed at leading order in an expansion in 1 /c can 
approximate quite well the numerical results. Moreover, 
a superposition of three QNMs matches the RD phase. 
In Sec. VI we apply the tools developed in the previous 
sections to understand, using analytic expressions, how 
the kick builds up during the inspiral, merger, and ring- 
down phases. Finally, Sec. VII contains a brief discussion 
of our mains results and future research directions. 


II. NUMERICAL SIMULATIONS 

In this section we introduce our definitions and no- 
tations, and review the main features of the numerical 
simulations. Throughout the paper, we adopt geometri- 
cal units with G — c— 1 (unless otherwise specified) and 
metric signature (— 1 , 1 , 1 ,!). 


A. Definitions and conventions 


Our complex null tetrad is defined using the time-like 
unit vector normal to a given hypersurface f, the radial 
unit vector r and an ingoing (£) and outgoing (n) null 
vectors as 
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We define the complex null vectors m and m* by 

m = -L((9 + i^), (3) 

m* = (4) 

with the standard spherical metric at infinity ds 2 = 
-dr 2 + dr 2 + r 2 (d6 2 +sin 2 Qdp 2 ). In terms of this tetrad, 
the complex Weyl scalar ^4 is given by 

^4 = Cabcd n a (m b )*n c (m d )* , (5) 

where C a bcd is the Weyl tensor and * denotes complex 

conjugation. 

To relate #4 to the gravitational waves (GWs), we note 
that in the transverse-traceless (TT) gauge, 
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Following usual convention, we take the /i+ and h x po- 
larizations of the GW to be given by 


h+ = \(h) J-hTT), (7) 

*x = (8) 


Since the Riemann and Weyl tensors coincide in vacuum 
regions of the spacetime: R a bcd = C a bcd , we find 


* 4 = -(ii + -i&x). (9) 


It is most convenient to deal with ^4 in terms of its 
harmonic decomposition. Given the definition of ^4 in 
Eq. (5) and the fact that m* carries a spin- weight of — 1, 
it is appropriate to decompose ^4 in terms of spin- weight 
—2 spherical harmonics - 2 YemiPiV) [47]. There is some 
freedom in the definition of the spin- weighted spherical 
harmonics. Here, we define them as a linear combination 
of the scalar spherical harmonics Y^m and Y^_i) m [48], 
that is 
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for i > 2 and |m| < l, and with the functional coefficients 
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Finally, in the far field (r » M) we decompose the di- 
mensionless Weyl scalar Mr $4 as 

MrVf 4 (t,r) = -2C em (t)~2Y im (9, p) . (12) 

tm 

where M = mi + m 2 is the total mass of the binary 
system with mi and m2 the masses of the individual BHs 
(see below for explanations), and r is the radial distance 
to the binary center of mass. 


The mass and spin parameters of the final BH are Mf 
and af. The values of Mf and af listed in Table I are 
computed from the loss of energy and angular momentum 
from the initial time to the end of the RD phase. They 
are compatible with the values obtained by extracting 
the fundamental QNM (see below Sec. IV B). All spins 
are orthogonal to the orbital plane, so = A y = 0 (the 
exception is a single run with planar spins discussed in 
Sec. VI D). 


B. Details of numerical simulations 

We set up the simulations by placing the BHs on an 
initial Cauchy surface using the Brandt-Briigmann pre- 
scription [49]; the Hamiltonian constraint is solved using 
the second-order accurate multigrid solver AMRMG [50]. 
We use the Bowen- York [51] framework to incorporate 
the BH spins and momenta, with the choice of initial 
tangential momentum informed by the quasi-circular PN 
approximation of Damour et al. [32], Eq.(5.3). 

The parameters for the runs considered in this paper 
are shown in Table I. We use the following notation. We 
indicate with EQ and NE the equal-mass and unequal- 
mass runs. The subscripts 0, +, — refer to zero spin, spin 
aligned, and spin anti-aligned with the orbital angular 
momentum, respectively. For the unequal-mass cases we 
use a superscript to indicate the mass ratio mi : m2. We 
denote by mi the BH horizon mass computed as 
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where Si = aimiSi = SiS\ is the spin angular mo- 
mentum of hole 1, ruir^i is its irreducible mass m* rrj i = 
a/Ai/ 167 t [52], and A\ is its apparent horizon area. Sim- 
ilar definitions hold for hole 2. The binary’s total mass 
is M = mi + m2, 5m = m\ — m2, the mass ratio 
is q = mi/ m2 < 1, and the symmetric mass ratio is 
77 = mim2/M. Following Kidder [35], we further define 
the spin vectors S = Si + S 2 , A = M(S 2 /m 2 - Si/mi), 
and ^ = S — {5m /M) A. 


The simulations were carried out using the moving 
puncture method [53, 54] in the finite-differencing code 
Hahndol [55], which solves the Einstein equations in a 
standard 3+1 BSSN conformal formulation. Dissipa- 
tion [56] terms (tapered to zero near the punctures) and 
constraint-damping [57] terms were added for robust sta- 
bility. We used the gauge condition recommended in 
Ref. [58] for moving punctures, fourth-order accurate 
mesh-adapted differencing [59] for the spatial derivatives, 
and a fourth-order accurate Runge-Kutta algorithm for 
the time-integration. The adaptive mesh refinement and 
most of the parallelization was handled by the software 
package Paramesh [60], with fifth-order accurate interpo- 
lation between mesh refinement regions. 

The grid spacing in the finest refinement region around 
each BH is hf = 3M/160. We extract data for the ra- 
diation at a radius r ex t — 45 M. The wave extraction 
was performed by 4th order interpolation to a sphere fol- 
lowed by angular integration with a Newton-Cotes for- 
mula. We have found satisfactory convergence of the re- 
sults. For example, for the 2:1 mass ratio run, for which 
a higher resolution of hf = 1M/64 was run in addition to 
hf — 3M/160, the rates of convergence of the Hamilto- 
nian and momentum constraints are comparable to that 
found in our equal mass runs reported in [61], and the ra- 
diated momenta from the two resolutions agree to within 
2%. This was also true for a 3:2 mass ratio test case with 
spin (NE++ below). 
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TABLE I: Parameters of the numerical simulations (see Sec. II B for explanations). 


Run 

mi m 2 5m q 

ai/mi a2/m2 

A z 

r 

Sl 3 

M f 

af/Aff 

EQ + _ 

0.503 0.503 0.0 1.0 

0.198 

-0.198 

-0.2 

0.0 

0.075 

0.9668 

0.697 

NEgi? 

0.401 0.593 -0.192 0.677 

0.0 

0.0 

0.0 

0.0 

0.0 

0.9599 

0.675 

NEob 2 

0.333 0.667 -0.333 0.500 

0.0 

0.0 

0.0 

0.0 

0.0 

0.9662 

0.633 

NEji, 4 

0.2 0.8 -0.6 0.250 

0.0 

0.0 

0.0 

0.0 

0.0 

0.961 

0.423 

NE+i 

0.399 0.610 -0.210 0.655 

0.201 

-0.194 

-0.2 

0.002 

0.072 

0.9707 

0.640 

NEi ; ® 

0.399 0.610 -0.212 0.653 

-0.201 

0.193 

0.2 

-0.002 

-0.072 0.9667 

0.704 


III. MULTIPOLAR FORMALISM 

In this Section we review the most relevant results from 
Thorne [41], showing how a multipole decomposition of 
the gravitational radiation field can be used to calculate 
the energy, angular momentum, and linear momentum 
fluxes from a BH binary system. One of the fundamental 
assumptions we employ here is that it is possible to relate 
in a simple way the multipole moments of the source and 
the radiation field, so in much of the discussion below we 
will use these two descriptions interchangebly. 

A. Linear momentum flux 

In the literature [7-11, 15] it is common to compute 
the linear momentum flux, and then the recoil, using the 
following formula 



where r is the extraction radius and the antiderivative of 
TU is used because the linear momentum flux scales as the 
square of the first derivative of the waveform, whereas ^4 
is proportional to the second derivative of the waveform 
[see Eq. (9) above]. To study how the different multipole 
moments contribute to the recoil, we could plug Eq. (12) 
into Eq. (14), as done, e.g., in Ref. [10]. Here, we pre- 
fer to use the expression of the linear momentum flux 
given in terms of the symmetric and trace- free (STF) ra- 
diative mass and current multipole moments, as done in 
Refs. [41-45]. 


Starting from Eq. (4.20’) in Ref. [41], we can write the 
linear momentum flux as 
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where 1^ (S A t ) are the STF /-dimensional radiative mass 
(current) multipole moments and left-hand superscripts 
represent time derivatives. From these tensors, we can 
construct the radiative multipole moments X lm and S lm 
according to the normalization given by Eq. (4.7) of 
Ref. [41]: 


X lm 
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Im 
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where are /-dimensional STF tensors that are closely 
related to the usual scalar spherical harmonics by 

Y lm (e,<p)= (17) 

with ri 1 = (sin 9 cos <p, sin 9 sin ip, cos 9 ). Note that the ra- 
diative moments X lm and S lm are complex scalars and 
have no explicit spatial dependence. To simplify the no- 
tation below, we incorporate the (/ + 1) time derivatives 
into the radiative multipole moments, and define 

jlm __ glm __ ( 1 + 1 ) ^ rn ( 18 ) 

By combining Eqs. (15), (16), and (18), we find that at 
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leading order (in a 1/e expansion) the linear momentum flux is given by 


Fj 0) + iF W = — — 
x v 336tt L 
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Note that Eq. (19) coincides with Eq. (9) in Ref. [38] 
when we express the radiative multipole moments in 
terms of the source moments [42-45] and reduce ourselves 
to a circular, non-spinning orbit in the (x, y) plane. In 


this case only the first three terms in Eq. (19) survive. 

The next highest order terms (1/c 2 ) are proportional 
to the mass octupole 7 3m , or current quadrupole S 2m : 


F^+iF^ = 


1 

672tt 


| -7iV^S 32 I 33 * - 14a/6/ 33 / 44 * + 4V21S 20 S 31 * - 4 a/355 21 5 32 * - 2a/2105 22 5 33 *+ 
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Note that all of the terms in Eqs. (19) and (21) contain 
products of multipoles with m' = m ± 1, while the terms 
in Eqs. (20) and (22) have m! — m, as with familiar 
quantum-mechanical operators that involve similar Xi- 
weighted integrations over the sphere. Also note that for 
mass-mass and current- current terms, V = l ± 1. 

The above formulae (19)-(22) are valid for completely 
general orbits, including spin terms and even for binary 
systems precessing out of the plane. However, we can 
simplify them significantly by rotating into the frame 
where the instantaneous orbital angular momentum is 
along the 2 -axis. Furthermore, by relating the radiative 
multipole moments to the source multipole moments [42- 


45] and assuming that terms proportional to R (R be- 
ing the binary radial separation) are negligible, we find 
J 20 = I 30 = S' 30 = J 32 = 7 40 = 7 41 = 7 43 = 0. In the ap- 
proximation of R — 0, the inclusion of R ^ 0 terms adds 
no new multipole modes. In the case of non-spinning 
BHs, the formulae (19)-(22) can be additionally simpli- 
fied by setting S' 20 = 7 21 = S 22 = S 31 = S 33 = 0. Quite 
interestingly, we obtain that the latter conditions are also 
valid in the special case of non- precessing BHs where the 
spins are aligned or anti- aligned with the orbital angular 
momentum. Since these are the cases we consider in this 
paper, we shall refer often to the following approximate 
formula for the linear momentum flux: 


F x + iFy 


1 r 
672t r . 


— 28iS 21 7 22 * 


- 2v^l07 22 7 33 * - 14\/67 33 7 44 * + 2\/l47 31 7 22 * - 7iv / 65 32 7 33 *] , 


F z = 0. (23) 


As we will see below in Sec. IV A, the linear momentum flux contributions from J 31 / 22 * as well as other higher-/ 


modes are typically smaller by at least an order of mag- 
nitude. When integrating Eq. (23) to get the recoil ve- 
locity, we also find that (due in large part to the rel- 
ative phases between the modes) the contribution from 
5 32 /33* ra ther minimal. Thus for most of the analysis 
that follows, we will focus solely on the first three terms 
ofEq. (23). 

In the following, sometime we shall use 

F = {F x ,F y ,F z }, F = |j. (24) 

All the (non-precessing) numerical simulations we will 
analyze have F z = 0, so we can introduce a complex 
scalar flux 


F = F x + iF y . (25) 

Since what we extract from the numerical simulations 
are the modes _ 2 Ce m computed over the sphere surround- 
ing the binary, we need to relate the _2 Ce m to the radia- 
tive mass and current multipole moments defined above. 
From Eq.(4.3) of [41], 

hmm = X( (i) X' m T a f Mm m a rn i> + ^S lm T^ Jm m a m b ) 

lm 

(26) 

where h mrn = h a bm a m b and h a b is the metric perturba- 
tion g a b — rjab in the transverse traceless gauge, which 
satisfies Eqs. (6a-6b), and T^ 2,lm and T^ 2,lm are the 
“pure-spin” harmonics of Thorne. From Appendix A of 
[62] 


T ab' lm = ^(-2 Y lm m a m b + 2 Y lm fh a m b ) (27) 

T ab’ lm = ^(-2Y lm m a m b - 2 Y lm m a fh b ). (28) 

Substituting Eqs. (27-28) into Eq. (26) and recalling that 
m a m a = 0 gives 

hmm = i T( (l) T lm + i (l) S lm ) +2 Y lm (29) 

^ 2 lm 

Now taking the complex conjugate and using the fact 
that +2 y* lrn — (— l) m - 2 Y l ~ m [note there is a typo in 
Eq. (3.1) of Ref. [47]] we obtain 

hmm = 4= y'(-l) m ( ( °^ m - i {l) S* lm ) _ 2 F' -Tn (30) 

^ 2 lm 

= 4= £(-ir( (,)r '" m - i {l) S* l ~ m ) - 2 Y\ 51) 

V lm 

Using the tetrad choice of Eqs. (1-9), d 2 h m m = h + — 
ih x = —^4 which decomposed into spin -2 weighted har- 
monics gives 


d 2 t hmm = - X) -2 C lm -2 Y lm , (32) 

lm 


we then see term-by-term that 

(_l)m(((+2)j*(-m _ i (!+2) <s *;-m) = _^2_ 2 Q m (33) 

Recall that (_1)™X*-™* = X lm and (-l) m 5 ! - m * = S lm , 
which allows us to write 

(!+ 2)l' m = ~[- 2 C tm + ,(34a) 

(i+ 2 ) s im = (-1 ) m _ 2 c;_ m ] .(34b) 

Equations (19)-(23) are expressed in terms of I lm = 
(i+i)jrim an< ^ gim = (i+i)gim w hi c h can be computed by 

integrating Eqs. (34a), (34b). To avoid the complication 
of an undetermined constant of integration, we typically 
integrate - 2 GmM backwards in time, since in the nu- 
merical data (and what we expect happens in reality) all 
the moments go to zero exponentially after the merger. 
At early times on the other hand, most of the modes are 
significantly non-zero and also include a large amount 
numerical noise due to the initial conditions. 


B. Energy and angular momentum flux 


Unlike the equations for the linear momentum flux, 
which all involve “beating” between pairs of different 
modes, the energy and angular- momentum flux expres- 
sions involve terms of the form \I lm \ 2 . As we will see 
below, for the comparable-mass binary systems that we 
analyse (mi : m 2 = 1 : 1, 2 : 3, 1 : 2), the amplitude 
of the mass quadrupole moment 7 22 is roughly an order 
of magnitude larger than the next largest mode. Thus 
it almost completely dominates the energy and angular 
momentum fluxes, and we can write [see Eq. (4.16) in 
Ref. [41]] 


dE 

dt 


32tt 


jr x (i/' m i 2 + is' m i 2 ) 


1—2 m——l 


167T 


1 1 


22 1 2 


(35) 


For the numerical simulations considered in this paper, 
the only non-zero modes have l + m even for I lm and 
l + m odd for S lm , so we can neglect the (m, m ± 1) 
cross-terms in Eq. (4.23) of Ref. [41]. These cross-terms 
are responsible for angular momentum loss in the x — y 
plane, so it is reasonable that they must be zero for non- 
precessing planar orbits. In this case, where the angular 
momentum is in the z direction, we have 


dJz 

dt 


i 

327T 


00 1 

X y V +1 )l lm + 

1=2 m=—l 


( l ) glm* (l+l) glm 

J_ (2)j22* (3)j22 


(36) 


where in the second equation we have restored the ex- 
plicit time derivatives [see Eq. (18)]. 
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TABLE II: Energy and angular momentum radiated in each of the dominant multipole modes. In parentheses we show the 
amount radiated only after the peak of each mode. 


Run 

E22 

£21 


£32 


£33 


E44 


J 22 

J21 


J 32 


J 33 


J 44 


EQ + _ 

0.035 

2.2 • 

■10" 

-5 

1.6- 

10' 

-4 

3.8- 

10- 

-6 

3.3' 

' 10- 

-4 

0.22 

7.0 ■ 

10“ 

-5 

7.9' 

• 10- 

-4 

2.0 

• 10- 

-5 

1.9- 

10“ 

-3 


(0.014) 

( 1 . 7 . 

■10" 

- 5 ) 

(1.2- 

10' 

■ 4 ) 

(2.3. 

10- 

- 6 ) 

( 1 . 5 . 

' 10' 

- 4 ) 

(0.050) 

(4.6 • 

10- 

- 5 ) 

(2.0. 

' 10- 

- 5 ) 

(1.2 

• 10- 

- 5 ) 

(6.4 • 

10' 

- 4 ) 

NEqo 3 

0.031 

6.1 « 

• 10' 

-5 

9.0- 

10' 

-5 

5.6* 

• 10- 

-4 

2.9' 

' 10' 

-4 

0.22 

2.1 ■ 

10- 

-4 

3.9- 

• 10- 

-4 

3.1 

• 10- 

-3 

1.8- 

10' 

-3 


(0.011) 

(4.0' 

• 10" 

- s ) 

(6.6- 

10- 

- 5 ) 

(2.8- 

10" 

- 4 ) 

(1.0 . 

• 10- 

- 4 ) 

(0.045) 

(9.8 ■ 

10“ 

■ 5 ) 

(2.5. 

10- 

■ 4 ) 

( 1.1 

• 10- 

- 3 ) 

(4.6 • 

10" 

• 4 ) 

NEjx) 2 

0.025 

1.4- 

' 10' 

-4 

4.7* 

10' 

-5 

1.2- 

10" 

-3 

2.7- 

10- 

-4 

0.18 

4.8- 

10- 

-4 

2.4- 

10- 

-4 

6.9 

• 10- 

-3 

1.7- 

10- 

-3 


Fo 

Lj 

0 

1 

CO 

(9.4- 

10' 

- 5 ) 

(3.0- 

10" 

• 5 ) 

(5.8- 

10- 

- 4 ) 

(7.3 • 

• 10- 

■ 5 ) 

(0.037) 

(2.4 • 

10" 

- 4 ) 

(1.3- 

10' 

- 4 ) 

(2.3 

• 10- 

- 3 ) 

(3.0- 

10- 

- 4 ) 

NEqo 4 

1 

0.012 

2.1- 

10" 

-4 

2.7- 

10- 

-5 

1.6- 

10- 

-3 

3.3 

' 10- 

-4 

0.12 

8.0- 

10- 

-4 

1.6- 

10- 

-4 

0, 

.011 


2.4- 

10- 

-3 


(3.5 • 1CT 3 ) 

(1.4- 

10- 

- 4 ) 

(9.3- 

10- 

- 6 ) 

(6.6- 

10“ 

- 4 ) 

(1.2. 

10' 

- 4 ) 

(0.016) 

(3.8 ■ 

10- 

- 4 ) 

(2.7 ■ 

10- 

■ 5 ) 

(2.9 

• 10" 

- 3 ) 

(4.8- 

10- 

• 4 ) 

NE+i 

0.029 

1.6' 

10- 

-4 

9.3- 

10- 

-5 

5.2- 

10" 

-4 

2.6. 

10" 

-4 

0.20 

5.4 ■ 

10" 

-4 

2.1 • 

10" 

-4 

2.9 

• 10" 

-3 

1.6- 

10- 

-3 


(0.010) 

(10 • 

10' 

- 4 ) 

(6.7- 

10- 

■ 5 ) 

(2.5 • 

10- 

- 4 ) 

(8.2 - 

•10" 

• 5 ) 

(0.031) 

(2.9 • 

10" 

■ 4 ) 

(5.3- 

10" 

- 4 ) 

(9.8 

• 10" 

- 4 ) 

(3.3- 

10" 

- 4 ) 

NEl% 

0.033 

1.4 • 

10“ 

-5 

1.1 • 

10- 

-4 

7.1 ■ 

10" 

-4 

2.9- 

10" 

-4 

0.23 

5.0 • 

10" 

-5 

4.4- 

10" 

-4 

3.9 

• 10" 

-3 

1.8- 

10- 

-3 


(0.011) 

(8.7- 

10' 

-6) 

(7.8 • 

10" 

-5) 

(3.4- 

10' 

- 4 ) 

(9.2- 

10- 

■5) 

(0.044) 

(2.1- 

10“ 

■ 5 ) 

(3.1- 

10" 

• 4 ) 

(1.3 

• 10" 

- 3 ) 

(3.7- 

10- 

■ 4 ) 


Integrating Eqs. (35) and (36) term-by-term, we can 
calculate how much energy and angular momentum are 
radiated in each of the dominant modes, similar to the 
approach of [68]. These results are shown in Table II, 
along with the contributions from just the inspiral phase 
(t < tpeak) and just the RD phase ( t > t pe a k). We will 
see below in Section V that these various energy con- 
tributions agree closely with the Newtonian predictions 
for the relative mass-scalings. For example, the energy 
E 22 in the inspiral phase should scale as 77, while the 
RD contribution should scale like rj 2 . It is important to 
note that the different moments have different scalings: 
Ess ~ 5m 2 , while the 7 44 contribution has a much weaker 
dependence on mass ratio: £44 77 2 (1 - 3t7) 2 - 

In the limit of very large initial separation (small ini- 
tial frequency), each of the Ei m and Jj m should converge 
to a finite value, with the notable exception of J 22 . It is 
well-know that the angular momentum of a binary sys- 
tem scales as I? 1 / 2 , and is thus unbound in the limit of 
R 00, but it is interesting to see that the higher-order 
contributions to the angular momentum all converge at 
large R . This can be understood directly from Eq. (36) 
in the Keplerian limit of R = o;“ 2 / 3 . At leading order, 
radiation reaction follows the relation dt ~ u~ ll ^du so 
the angular momentum in the inspiral is 

J 22 = ± f t0 dt^l 22 * (3) J 22 

y t=B _ 00 

rv 0 

~ / uW^o;- 11 / 3 - 00. (37) 

J Lj = 0 


For all the other energy and angular momentum 
modes, the fluxes from Eqs. (35), (36) scale as a; 10 / 3 or 
higher powers, and thus converge when integrated over 
u- n /*du. 


IV. MULTIPOLE ANALYSIS OF THE 
NUMERICAL SIMULATIONS 

In this section we want to investigate how the different 
multipole moments contribute to the inspiral phase and 
the RD phase of unequal-mass binary systems. 


A. Inspiral phase 


As it can be derived in PN theory [27] and it has been 
confirmed numerically in Refs. [53, 54], the l = 2,m = 2 
mode in Eq. (12) is circularly polarized during the in- 
spiral phase of a non-spinning equal-mass binary system. 
Because of that, Ref. [66] defined the (dominant) orbital 
angular frequency as 

<38 » 

Here, we extend (38) by defining several (dominant) 
angular-orbital frequencies, each of them being related 
to a specific multipole moment, I lm or S lrn , as 


UJ 


Ilm 

D 



U 


Sim 

D — 




• (39) 


F 




FIG. 1: Dominant orbital angular- frequency obtained from the individual radiative multipole moments, as determined by 
Eq. (39). The different frequencies with l = m agree well during the inspiral phase but start decoupling during the transition 
to RD. The frequency with l = 2,m = 1 decouples from the others at earlier time and reaches a much higher plateau. The 
left panel refers to the NEqo 3 run and the right panel to the NEjo 2 run. We denote with t pea k the time at which I 22 reaches its 
maximum. 




(t-W)/M HeJ/M 

FIG. 2: Amplitudes of the dominant radiative multipole moments. On the left panel we show the modes for the NEqo 3 run, 
while on the right panel the modes for the NEqo 2 run. The leading-order mass quadrupole 1 22 is about an order of magnitude 
stronger than any other mode. The oscillating behavior of the S 32 moment is likely due to mode coupling. We denote with 
£ P eak the time at which I 22 reaches its maximum. 


We plot these frequencies in Fig. 1 for the dominant 
multipole moments, i.e., 7 22 , S' 21 , I 33 , 7 44 , and S 32 . The 
amplitudes of the I 31 and I 42 modes are too weak and 
dominated by noise to extract a dominant frequency. In 
this figure, as well as most shown in the rest of the paper, 
we plot the time variable with respect to t pe a k, the point 
at which \I 22 \ reaches its peak, closely corresponding to 
the peak in GW energy emission [see Eq. (35)]. We notice 
that the frequencies associated to the modes with l = m 
agree quite well during the inspiral, but for the frequency 
of the S ' 21 mode decouples from the others approximately 
50 M before the peak in the 1 22 mode. As we shall see in 
Sec. VI, this is due to the fact that, during the ringdown 
phase, the dominant angular frequency associated to the 
S ' 21 mode is almost twice as large as those of the other 


leading modes [63-65]. This decoupling plays a major 
role in determining the shape of the kick and anti-kick 
(see Sec. VI below). 

In Fig. 2 we show the amplitudes of the multipole mo- 
ments in Eq. (23). The’ left panel in Fig. 2 refers to the 
NE 2 5 3 run, while the right panel to the NEJ 5 2 run. The 
mass-quadrupole moment I 22 clearly dominates in both 
cases, while the I 31 and 7 42 modes are so weak as to 
be almost completely overwhelmed by numerical noise. 
In addition to having dissimilar amplitudes, the differ- 
ent moments also peak at slightly different times, which 
may be related to the fact that RD modes are excited 
at different times. In particular, the modes mentioned 
above with l ^ m tend to peak later in time, perhaps 
due to a longer transition to the higher QNM frequency. 
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(t-W)/M 



HeakVM 


FIG. 3: Linear momentum flux of the strongest radiative multipole moments, i.e., the ones in Eq. (23). On the left panel we 
show the modes for the NEqo 3 run, while on the right panel the modes for the NEqo 2 run. We denote with t pea k the time at 
which / 22 reaches its maximum. 


As we shall see in Sec. V, as the mass ratio becomes more 
extreme, the higher-order modes increase in relative am- 
plitude, with I 33 and S 21 both proportional to 5m. I 44 
and 5 32 , however, scale as 77(1 — 377 ), so they increase 
only slightly in the range of masses considered here. 

Next, in Fig. 3, we show the amplitude of the linear mo- 
mentum flux from the mode-pairs included in Eq. (23). 
Again, the mass-quadrupole terms dominate, with signif- 
icantly smaller contributions from the 5 32 and 1 31 modes. 
However, note the appreciable flux amplitude from the 
j33j 44* t erm? w hich is formally a higher-order correction 
in a (1/c) expansion [37, 38]. From Fig. 3, we expect 
that the first three pair of modes in Eq. (23) should con- 
tribute most significantly to the recoil. Including the 
complex phase relations between the different modes, we 
find this result will be supported further by the analysis 
in Sec. VIA. 


B. Ringdown phase 

We now extract the QNMs, notably the fundamental 
and the first two overtones, present in the most significant 
multipole moments. We follow the procedure outlined in 
Ref. [66]. To avoid possible constant offsets introduced 
by integrating Eqs. (34a), (34b), we prefer to extract the 
QNMs directly from the _ 2 C f / m ’s instead of using Ii m 's 
or Sim's. Following the approach of [65], we define the 
complex frequencies cri mn : 

&lmn — & Imn ^ / ^Imni (40) 

and each RD mode is proportional to exp(— icri mn t). In 
this format, ui m n are the QNM oscillation frequencies 
(not to be confused with the dominant frequencies of 
Eq. (39)) and r/ m n are the mode decay times, all func- 
tions of the final black hole mass and spin. In this nota- 
tion, l and m are the same spherical wavenumbers used 


above, and n = 0 denotes the fundamental mode, with 
n = 1, 2, . . . corresponding to the higher harmonics. The 
fundamental QNM frequencies cq m 0 are shown in Ta- 
ble III for the NR runs listed above. All frequencies and 
decay times are measured in units of the final mass Mf. 


We present the RD analysis only for the NEqq 3 run, but 
the others are qualitatively very similar. By applying a 
4-parameter fit to the most dominant mode, i.e., -2C22, 
we obtain af/Mf = 0.669 and M/Mf = 0.965 together 
with the amplitude and phase of the fundamental QNM. 
Using the results of Ref. [65], we find Muj 220 = 0.538 and 
Af /7220 = 0.0847. We then do a 6-parameter fit and suc- 
cessfully extract also the first overtone, obtaining slightly 
different values for af/Mf = 0.661 and M/Mf = 0.958. 
We find it impossible to extract, with a 8-parameter fit, 
also the second overtone. By contrast if we keep af /Mf 
and M/Mf fixed and equal to the values obtained when 
extracting the fundamental QNM, we obtain that we can 
fit up to the second overtone. Moreover, quite interest- 
ingly, the fit provides a waveforms which compares very 
well with the NR waveforms up to ~ 5M before the peak 
of - 2 ^ 22 ? as can be seen in the upper left panel of Fig. 4. 

The other panels of Fig. 4 show results for the other 
relevant modes - 2 C 33 , - 2 C 44 and - 2 C 32 . As obtained in 
Ref. [66], we find a mode-mixing in - 2 ^ 32 , be., the QNM 
is a combination of l = 2 , 777 = 2 and l = 3, m = 2. By 
fitting the fundamental QNM we obtain af/Mf = 0.670 
and M/Mf = 0.972; a f /M f = 0.412 and M/Mf = 0.838; 
af/Mf = 0.628 and M/Mf = 0.972, for - 2633 , - 2 C 44 and 
- 2 C 32 , respectively. Those values for the final BH spin 
and mass are rather consistent, except for - 2 C 44 . This 
discrepancy might be due to numerical resolution effects, 
and will be the object of future investigations. 




TABLE III: Frequencies and decay times for the fundamental QNMs for each of the numerical simulations. 


TO 


Run 

af / Mf 

cr 210 

T210 

0220 

T220 

<7320 

T320 

C330 

T330 

C440 

T440 

EQ + _ 

0.697 

0.454 

12.2 

0.531 

12.4 

0.758 

11.9 

0.841 

12.0 

1.14 

11.8 

NEgb 3 

0.675 

0.450 

12.1 

0.521 

12.2 

0.749 

11.7 

0.827 

11.9 

1.12 

11.7 

NEJo 2 

0.633 

0.442 

11.9 

0.505 

12.1 

0.734 

11.6 

0.803 

11.7 

1.09 

11.5 

NEJi) 4 

0.423 

0.411 

11.5 

0.445 

11.5 

0.674 

11.1 

0.711 

11.1 

0.963 

10.9 

NE+i 

0.640 

0.443 

11.9 

0.507 

12.1 

0.736 

11.6 

0.806 

11.7 

1.09 

11.5 

ne! : + 

0.704 

0.456 

12.2 

0.533 

12.4 

0.760 

11.9 

0.845 

12.1 

1.14 

11.9 




(t-W)/M 




FIG. 4: Comparison of numerical and QNM waveforms for the NEqo 3 run. The dominant modes analyzed are -2C22 ( upper 
left), _ 2C33 ( upper right), -2C32 ( lower left), and -2C44 ( lower right). We denote with t pea k the time of the peak of I 22 . 


V. EFFECTIVE NEWTONIAN MODEL 


In an attempt to better understand the amplitudes and 
frequencies of the various modes during the inspiral and 
merger phases, we present here what we call the “effective 
Newtonian” (eN) model. It begins with calculating the 
leading-order Newtonian formulae for each multipole mo- 
ment of the source, as a function of the BH masses, binary 
separation, and orbital frequency. To extend these for- 
mulae through the end of the inspiral and into the merger 
phase, we introduce an effective radius coordinate that is 


a way of absorbing PN effects into the leading-order mul- 
tipole expressions. Each multipole moment is then indi- 
vidually matched to a linear superposition of ring down 
modes, as is done in the effective-one-body model. Taken 
together, this eN model provides an excellent framework 
within which we can understand the details of the linear 
momentum flux and integrated recoil velocity. 
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A. Newtonian Multipole Moments 

Working at leading Newtonian order for each mode, 
we equate the radiative multipole moments to the source 
multipole moments. Restricting ourselves to circular, 
planar orbits, and setting M = 1, we find that for non- 
spinning systems, the dominant modes are [42-45] 


r*21 

^nospin 

-i^\[^-SrnrjR 3 uj 4 e , 
3 V 5 

(41a) 

7-22 __ 
-‘nospin 

q 6 y / | r?i? 2 w 3 e -2# j 

(41b) 

I 31 . = 

nospin 


(41c) 

c32 _ 

° nospin 

- T 'Jy f (! - 3 i) " 5 e_2 ‘** 

(41d) 

7-33 

-‘nospin 

54 Sm rj R 3 uj 4 e~ 3 ^, 

(41e) 

7 42 . = 

-‘nospin 

^i\/27T7] (1 — 3 r}) R 4 uj 3 e” 2 ^, 
63 

(41f) 

7-44 __ 

1 nospin 

-xVr^ 1-3 ^ r4ljS e ~ u<t> 

,(41g) 


where R is the radial separation and uj = <j) is the binary 
orbital frequency. Considering only the mass quadrupole 
terms in the linear momentum flux, we obtain the well- 
known result valid at Newtonian order [38]: 

F (0) = -i ^<5m rj 2 R 5 uj 7 e**. (42) 

Including the next- highest order moments in Eq. (21), we 
get 

F (1) = ( 1 ~ 3?7)-R 7 w 9 e # . (43) 

While there may also be next-to-leading order contribu- 
tions from a PN expansion of the multipole moments in- 
cluded in Eq. (19) that would show up in Eq. (43), we can 
effectively absorb those corrections into the R variable, 
as will be descibed below. 

Combining Eqs. (42) and (43) we find the linear mo- 
mentum flux scales like 


|F(°) + F ( i)| oc 6m rj 2 


1 + Si< 1 - 3 ^ V 


^6mr] 2 (l - 0 . 977 ), 


(44) 


which is remarkably similar to the result found in [9]. 
Here we have used R 2 uj 2 « 0.23 — 0.25 at the peak of the 
energy flux, which seems to be quite standard for a range 
of mass ratios. However, the extremely close agreement 
with Ref. [9] is probably to some degree a coincidence, 
since this simple Newtonian formula does not include any 
details of the phase relations between different modes. 
As we show below in Sec. VI B, in the transition from 
inspiral to ringdown, these phases play a major role in 
determining the amplitude of the un-kick, and thus the 
final recoil velocity. We find that more extreme-mass 
ratio BH binaries have a relatively smaller un-kick, which 
should also play an important role in the scaling relation 
of [9]. 


If we compute the above multipole moments (41a)- 
(41g) using uj as given by Eq. (39) and R as obtained 
from the puncture trajectories, we do not find a very good 
agreement with the numerical results. This is not sur- 
prising since there is no reason to believe that the New- 
tonian approximation should work well all along the in- 
spiral phase. We should expect that higher-order PN cor- 
rections become important as we approach the merger. 
Since our scope is limited to a diagnostic of the NR re- 
sults, and not to a precise comparison with PN calcula- 
tions, instead of including PN corrections in Eqs. (41a)- 
(41g) and (42), we investigate whether by properly scal- 
ing the Newtonian expressions we can get a better agree- 
ment until the merger. We can also think of this normal- 
ization as a way of resumming the PN expansion. 


Quite interestingly, if we compute the amplitudes | I lm \ 
or |5 Zm | from the numerical data, and the angular fre- 
quency uj from Eq. (39), we find that the radii Ri m which 
appear in the RHS of Eqs. (41a)-(41g) are rather inde- 
pendent of the multipole moments l and m, as Fig. 5 
shows. We denote the radii Rim computed numerically 
as effective radii Rf^. The close agreement between the 
frequencies (see Fig. 1) and effective radii for each mode 
suggests we can use the Newtonian expressions and a 
single i? eff (£), e.g., the f° r modes with a high 

degree of accuracy for the entire inspiral phase and even 
during the transition to merger. 


For comparison we show in Fig. 5 also the radius 
from the puncture trajectory (dot-dashed curves) and 
the radius computed using the Arnowitt-Deser-Misner 
transverse-traceless gauge (dashed curves) , given as func- 
tion of frequency through 3PN order by [67] 


Radm = w 2 ^ 3 


1+uj 


2/3 


(- i+ d 


■ UJ 


.4/3 


+ o^ + ^r + w -7 - 


1 

’4 


1625 167 o 3 o 2 , 

77 -\ riTt 77 H m 

144 ' 192 1 2 81 ' 


( 45 ) 


TZ 




FIG. 5: Effective radius for different modes, derived from Eqs. (39), (41a)-(41g). The close agreement for the suggests 
we can use a single effective radius R eff (£) for the Newtonian expressions. We believe that the large oscillations in are due 
to initial eccentricity at early times. Also plotted is the ADM radius (dashed curves) derived from the orbital frequency via 
Eq. (VA) and the coordinate separation of the BH punctures (dot-dashed curves). The results correspond to the NEqo 3 (left 
panel) and NEqo 2 (right panel) runs. We denote with t pea k the time at which / 22 reaches its maximum. 


In the next section, we shall investigate how this simple 
effective Newtonian model can be combined with a su- 
perposition of QNMs, as described in Sec. IV B, giving a 
good representation of the NR results. 


B. Matching to ringdown 


We now match the inspiral and RD waveforms in a 
mode- by-mode fashion following the philosophy of the 
EOB approach [30]. A similar attempt was followed in 
Ref. [38] , where for simplicity the authors performed the 
matching to the Schwarzschild QNM frequencies, while 
we use the Kerr QNM frequencies and match to the fun- 
damental QNM frequency and the first two overtones, as 
done in Ref. [66]. We obtain the QNM frequencies and 
decay times from Ref. [65] . For the fundamental and two 
overtone QNMs, we can match a given multipole mode 
by equating it and two time derivatives to a linear com- 
bination of QNMs. 

We write 


I lm (t) = A(t) e ^(0 — Ai mn e i(Tlrnn * match \ (45) 

n=0 


where the complex QNM frequencies are known functions 
of the final BH mass and spin, and we must solve for the 


complex amplitudes Ai mn . Matching three QNMs we get 

2 

I (^match) = ^ Amni 

n — 0 


(^match) — ^ ^ ^ GlmnAmni 


n —0 

2 


j2 

(^match) = ^ ^ ®lmn Amm 

n=0 


(47) 


or as a simple matrix equation 


1 1 1 ^ 

i&lmO ^Iml 
2 2 2 
^ImO ®lm\ ®lm 2 J 


( AmO \ 


I Ami I 
\ Am 2 ) 


( I lm \ 

jlm 

^ jlm J 


(48) 

In Fig. 6, we compare the NR modes to the modes 
obtained by the effective Newtonian model defined in 
Sec. V A until £ m atch and by the superposition of three 
QNMs for t > t match* During the inspiral, the different 
moments are calculated according to Eqs. (41a)-(41g), 
using a single R e ff and wn determined from the / 22 mode, 
with the exception of the S 21 mode, where we instead 
use the higher frequency o;^ 21 . The QNM frequencies 
and decay time are obtained from the mass and spin of 
the final BH (see Tables I, III), so they are not com- 
puted a priori . We treat t matc h as a free parameter: if 
we stop the inspiral too early, the effective Newtonian 
mode amplitudes are still growing, so the sudden tran- 
sition to decaying RD modes prematurely reduces them. 
On the other hand, if the inspiral is continued too long, 
we tend to lose the important phase shifts between the 
modes that only begin during the transition to RD. This 
is particularly evident in the 7 44 mode, which undergoes 
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FIG. 6: Comparison of the effective Newtonian and NR radiative modes during inspiral, merger and RD phases. The data 
refer to the NEjo 2 run. We denote with £ pea k the time at which I 22 reaches its maximum. 





an unexplained 180° phase-shift around the transition to 
RD, and also decays at somewhat different rate than is 
predicted from QNM theory (see above, Sec. IV B). Mo- 
tivated by the results of Sec. IV B, notably by the fact 
that a superposition of three QNMs can fit very well the 
NR waveforms starting from the peak of the energy flux, 
we choose as best-matching point the peak of the energy 
flux. 

Having shown a reasonably close match for each of 
the radiative multipoles between the effective Newtonian 
model and the numerical data, it stands to reason that 
the total recoil calculated with this model should agree, 
as well. This is shown in Fig. 7, where we have also 
varied the matching point around t pea k- The significant 
range in final recoil velocities is due to the phase shifting 
of the different modes around t — £ pea k> most notably 
that of the I 44 mode described above. In Section VI B 
below, we will explore this phasing in greater detail and 
show how it affects the overall kick. At this point, we 
unfortunately do not have a clear understanding of the 
underlying cause of the de-phasing, but it may well be 
related to the slightly different times of transition from 
inspiral to ringdown for the different modes. Preliminary 


results also suggest that this de-phasing effect is reduced 
in more extreme-mass-ratio systems. 

VI. ANATOMY OF THE KICK 

In the above Sections, we have laid the groundwork 
for a multipolar analysis of the gravitational recoil, de- 
scribing the momentum flux as a combination of radiative 
multipole modes. Along with the psuedo-analytic models 
for the inspiral and ringdown phases, we can now give a 
detailed description of the “anatomy” of the kick, namely 
the way the different modes combine to produce a peak 
recoil velocity, followed by a characteristic un-kick and 
then asymptotic approach to the final value of the BH 
recoil. 


A. Contribution from different moments 

In Sec. Ill A, we showed how the individual multi- 
pole moments contribute to the linear momentum flux 
through the intregral of the #4 scalar [Eqs. (12), (14)]. 
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FIG. 7: Comparison of the effective Newtonian model and NR predictions for the recoil velocity for a range of inspiral- RD 
matching points. We denote with £ pea k the time at which 1 22 reaches its maximum. The data refer to the NEqo 3 run. 


Here, we want to determine exactly which modes we need 
to include in the multipole expansion Eq. (15) to get a 
good representation of the full recoil, and which are the 
pairs of modes in Eq. (23) that contribute most. 

By including only a select choice of terms in the $4 ex- 
pansion Eq. (12), we can calculate the linear momentum 
flux by direct integration of Eq. (14) and compare it with 
the predictions of Eqs. (19)-(23), in each case including 
only the appropriate moments. This is a good way of 
double-checking those lengthy equations term-by-term, 
and in practice we find excellent agreement, limited only 
by the numerical accuracy of the simulations. Similarly, 
we can use this method of truncated expansion to deter- 
mine which modes are necessary for calculating the recoil 
up to a given accuracy. The results of using higher and 
higher order multipolar moments are shown in Figs. 8, 9 
for the NEq 5 3 and NEJ5 2 runs. 

In the left panel of Figs. 8, 9 we show with a solid curve 
the integrated momentum kick from Eq. (14), with a 
dashed curve the contribution from terms up to l = 4, i.e., 
the ones obtained by Eq. (19) and (21), and with a dot- 
ted curve the contribution from the leading three terms 
in Eq. (23). We conclude that the linear momentum flux 
is dominated by the / 33 / 22 * ? /33 j 44*^ anc i 521/22* terms, 
which combine to produce the primary kick and anti-kick 
agreeing with the exact result within <10% throughout 
the entire merger. Note that the flux from the 5 32 / 33 * 
term, while not insignificant in Fig. 3, contributes almost 
nothing to the net recoil velocity. This is largely due to 
phase relations between the various modes during the 
transition from inspiral to ringdown, described below in 
Sec. VI B. 

In the right panel of Figs. 8, 9 we show the difference 
between the calculation obtained including terms up to 
l = 3, 4, 5, 6, and the exact result. It seems clear that we 
need modes up to and including l — 4 to get an accurate 
estimate of the recoil velocity. For more extreme mass 


ratios, higher-order moments become relatively more im- 
portant, but remain strongly sub-dominant to the l < 4 
modes [68]. 

To understand more clearly the relative contributions 
of the different modes to the total recoil, we will in- 
clude analysis of a few more simulations including non- 
precessing spins. As mentioned above in Sec. Ill A, non- 
precessing spins do not introduce any additional mo- 
ments compared to the non-spinning simulations, but 


simply modify the relative amplitudes of the different 
modes in Eq. (23) by adding the spin terms. Thus, once 
we determine how the spins modify the individual modes, 
we can use the same analysis for the spinning and non- 
spinning cases. 

Again setting M = 1 and relating the radiative multi- 
pole moments to the source moments, we get the leading 
order spin-orbit modifications to Eqs. (41a)-(41g) [see 
Eqs. (3. 14), (3.20) in Ref. [35] and Eq. (5.5) in Ref. [69]]: 

q21 

^SO — 

- 4 i 7 ^t)Su 3 e-‘*A», 

(49a) 

t 22 

1 SO ~ 

ji 

(49b) 

q 32 

^so — 


(49c) 

/-31 _ 

i SO — 


(49d) 

7-33 _ 
i SO — 


(49e) 

where we have introduced the spin vectors 


S31 = 

y<5mS + 1(11 - 39?7)A, 

(50a) 

£33 = 

|jmS + 5(i_ 577 ) a. 

(50b) 
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FIG. 8: In the left panel we show the net recoil kick, integrated from the linear momentum flux via Eq. (14) (solid curve), from 
all modes with l < 4 (dashed curve) and also limiting the modal composition of 4>4 to just the three dominant mode pairs in 
Eq. (23) (dotted curve). In the right panel we show the difference between the exact result and the #4 expansion Eq. (12), 
limited to l < 3,4, 5,6. The data refer to the NE§o run. We denote with t peak the time at which I 22 reaches its maximum. 
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FIG. 9: Same as Fig. 8, but for the NEjo 2 run. 


In all of the simulations considered here, where the di- Then Eqs. (41a)-(41g), and (49) give the linear momen- 
mensionless spins are equal (ai/rai = ^ 2 /^ 2 ) and point turn flux F = F x + iF y during the inspiral for each of the 

in opposite directions, \ z — 0, so we are left only with first three dominant terms in Eq. (23): 

the modification of S 21 due to A* and I 33 due to £§ 3 . 


77121,22 _ 

^insp 

iV 2 R 3 w 6 (2 5m R 2 lv + 3A Z ) e**, 

(51a) 

1722,33 _ 

"Mnsp 

— R 5 u 1 (5m + uj E 33 ) , 

(51b) 

1733,44 __ 
"Mnsp 

S?A 

— —irj 2 (1 — 377 ) R 7 co 9 ( 8m + uj £ 33 ) e 

(51c) 


While these flux formulae contain terms of various orders 
in a;, the effective Newtonian scaling of R should ensure 
that we are including all relevant PN terms. For con- 


sistency, we find that one must be careful to distinguish 
between and cj ^ 21 in Eq. (51a). 

The amplitudes of these fluxes are plotted in Fig. 10 
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FIG. 10: Relative amplitudes of the dominant multipole mode-pairs in the linear momentum flux Eq. (23). Also shown in the 
dashed curves are the eN model predictions for the flux amplitudes. We denote with t pea k the time at which 1 22 reaches its 
maximum. 


for the four runs NE^, NE+ 3 , NE 2 q 3 , and EQ + _. As 
seen in Table I, the NE 2:3 run has A z = 0.2 M 2 , while 
the NE^i run has A 2 = — 0.2M 2 , respectively, adding 
destructively and constructively with the 6 m term in 
Eq. (51a). This difference is clearly seen in the blue 
curves in Fig. 10. Also notable in these plots is the differ- 
ence in mode amplitudes from / 22 / 33 * ? due to a similar 
effect from the constructive/ destructive additions of 6m 
and £33 in Eq. (51b). As we see in Fig. 10, NE 2 5 3 appears 
to be the average of NE^ and NE^, while the flux from 
EQ + _ is strongly suppressed due to the 6m = 0 terms in 
Eq. (51c), leaving only the flux from A* = — 0.2M 2 and 
£| 3 = 0.075. 

In each panel of Fig. 10, we also plot with dashed lines 
the eN prediction for the various flux amplitudes. In al- 
most all cases, the eN flux is quite close to the NR results 
up to about 10M before t pe a k, when the effective radii of 
Fig. 5 begin to diverge. The amplitude differences near 
the peaks are comparable to those see in Fig. 6 for the 
NE 2 5 3 run. The notable exception is the F 21 ’ 22 flux from 
the NE^ and EQ + _ runs, where there may be addi- 
tional spin-spin terms that become important at these 
relatively low amplitudes. 


B. Transition to ringdown and the de-phasing of 
the multipole modes 


Since the flux vectors defined in the previous section 
will not generally be co-linear (and in the case of process- 
ing spins, not even co-planar), to understand the time 
evolution of the recoil velocity, we must first understand 
the phase relations between the different modes. From 
Eqs. (41a)-(41g) and (51a)-(51c), we see that during 
the inspiral phase, the individual moments and the re- 
sulting flux vectors evolve according to a single orbital 
phase </>, with F^p 2 pointing in the opposite direction of 

^insp 3 an d ^insp 4. However, as we can see from Fig. 1, 
as the binary evolves from inspiral to RD, the frequency 
(and thus phase) of the S 21 mode decouples from the 
other dominant modes. Upon closer inspection, we find 
that even the I 22 , 1 33 and I 44 modes deviate from each 
other enough to undergo a significant phase shift at the 
inspiral-RD transition. 


To quantify these effects, we define the phase differ- 
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ences: 


cos^ 2-3 

_ £21,22 . £22,33 
r insp * r insp > 

(52a) 

cos^ 2-4 

_ £21,22 . £33,44 
r insp * r insp > 

(52b) 

cos t/? 3-4 

_ £22,33 _ £33,44 
r insp ’ r insp ‘ 

(52c) 


Here we use the notation '0 m-m to describe the phase 
difference between two complex flux vectors, where m 
and ml correspond to the larger m - values of each mode 
pair that makes up the flux. These definitions are valid 


throughout the inspiral, merger, and ringdown phases. In 
the inspiral phase, we can see that for the unequal-mass 
runs, where 6m dominates Eqs. (51a)-(51c), we have 

COS ^irisp = cos tfnst = - 1 > COS ^insp = 1 • (53) 

For the EQ + _ run with 5m = 0, Eq. (51a) predicts that 
all phases have cos^insp = 1 during the inspiral. Dur- 
ing the RD phase, using Eq. (46), we can write the flux 
vectors and phase evolution in terms of the fundamental 
QNM frequencies ui m q: 


and 


t-,21,22 

_ ™21,22 
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r match 

t-i22,33 
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(54c) 
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cos^Id 4 


= cos[(o> 2 io - 2 w 220 + O >330 )(t - imatch) + ^matchl > 

= cos[(w 21 o — w 22 0 — O»330 + U)44o)(t — £ ma tch) + ^m^tchl ) 
= cos[(w 22 o — 2^330 + U>44o)(t ~ £ match) + Snatch] • 


(55a) 

(55b) 

(55c) 


Here $ ma tch is a phase offset determined at the transition 
from inspiral to ringdown. Quite interestingly, we find 
that for the range of BH spin parameters 0.5 < af/Mf < 
0.9, the linear combinations of frequencies in Eqs. (55a)- 
(55c) vary by less than 10%. Thus, if we compute the 
above expressions for the ui m0 corresponding to af/Mf = 
0.7, we have [65] 
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Even more intriguing, we find that for the unequal- 
mass simulations described above, the phase relations 
during the inspiral and RD are almost identical, regard- 
less of spin orientations. This can be seen clearly in 
Fig. 11, which plots cos'ip during inspiral, merger and 
RD for the different runs. The colinearity of the flux 
vectors is clear during the inspiral phase, and the sinu- 
soidal oscillations of the phases during RD agree well 
with the analytic predictions (plotted in dashed curves 
in Fig. 11). Since the analytic models are most reliable 
during the inspiral and RD phases (but have more dif- 
ficulty tracing the merger portion), we omit in Fig. 11 
the transition region of — 10M < (t - t pea k) < 10 M. The 
phase relations during inspiral are determined by Eq. (53) 


and during ringdown by Eqs. (56a)-(56c). Here we use 
a ^match about 20 M after t pea ^ to ensure that the phase 
relations are truly dominated by the fundamental QNMs, 
and thus Eqs. (56a)-(56c) are valid. Note that the phase 
differences for EQ + _ are particularly noisy since the am- 
plitude of the / 33 moment is zero to leading order, and 
it is thus more difficult to extract a clear phase for that 
mode. 

The feature that is most difficult to explain from an 
analytic model alone (and is thus omitted from the eN 
curves in Fig. 11) is the 180-degree jump in phase be- 
tween i^nsp 3 an d ^Insp 4 aroun d —20 M before the peak. 
This appears to be a general feature in all the unequal- 
mass runs, but preliminary results suggest that is is less 
significant (i.e., a smaller phase shift) for more extreme- 
mass ratio systems. We are not able to explain it with 
the additional RD overtone modes described in Sec. VB, 
but using slightly different RD matching points for the 
different multipoles may help explain the issue. 


C. The anti-kick 

These flux amplitudes and phase relations can now be 
used to understand the amplitude of the kick and anti- 
kick. Throughout the inspiral phase, the amplitude and 
rotational frequency of the flux vectors in Eqn. (51a) are 
monotonically increasing, giving the familiar outward- 
spiraling trajectory for the velocity vector. Then, in 
the RD phase [Eqs. (54), (55)], the frequency is nearly 
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FIG. 11: Phase differences between different mode-pair flux vectors, as defined by Eqs. (52a)-(52c). The data refer to the 
NE+ 3 (upper left panel), NE^ (upper right panel), NEqo 3 (lower left panel), and EQ+_ (lower right panel) runs. The dashed 
curves are the eN model predictions of Eqs. (53), (56). We denote with t pea k the time at which / 22 reaches its maximum. 


constant while the mode amplitudes decay exponentially, 
giving an inward-spiral that decays like a damped oscilla- 
tor around the final asymptotic recoil velocity. These tra- 
jectories in velocity space can be seen in Figs. 12, along 
with the instantaneous flux vectors from the competing 
mode-pairs. Clearly, even small changes in the mass ra- 
tios and spins orientations of the BHs can give a rather 

1 


diverse selection of velocity trajectories. 

To calculate the recoil velocity, we must integrate the 
linear momentum flux vectors in time. We can get a 
reasonable analytic approximation by using Eqs. (51a) 
and (54) for the inspiral and RD phases, respectively. In 
the adiabatic inspiral, the complex recoil velocity v = 
v x -f iv y can be written as 


/ t match 2 

F{t')dt' ~ -F^match), (57) 

-oo ^match 


while for the RD portion we have 
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(58) 


summing over each pair of modes (Zm, I'm'), Then the total velocity v in each of the dominant mode pairs is given 
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FIG. 12: The recoil velocity vector evolving in the v x -v y plane (black solid curve), along with the flux vectors due to the three 
mode pairs at each time interval (same color scheme as Fig. 10). The data refer to the NE+i (upper left panel), NEi : + (upper 
right panel), NEqo 3 (lower left panel), and EQ+_ (lower right panel) runs. We denote with the label peak the time at which I 22 
reaches its maximum. 
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The phase ^atch defined as the angle made between 
the flux vector F 21 ' 22 and the integrated velocity vector 
v at the beginning of the ringdown (with other phases 

^match’ '/'match defined analogously). Because of the 
anomolous phase shifts and departure from adiabatic- 
ity at the transition from inspiral to ringdown, these 
angles are difficult to predict with an independent an- 
alytic model, but can be calculated easily from plots like 
Fig. 12. 

In the case of the NE+i run, where the recoil is almost 


entirely dominated by the mode pair / 22 / 33 * ? the ratio 
of the final velocity to the peak velocity (i.e., the ampli- 
tude of the anti-kick) is given by the magnitude of the 
complex number |1 - ia; matc he i7r/2 / ( 0*220 - 0 * 330 )! ~ 0.5. 
For the EQ + _ run on the other hand, we see that there 
is NO anti-kick, which can be explained by the slowly ro- 
tating flux vector that does not spiral back inwards, but 
rather drifts off slowly towards infinity during the RD. 
The difference between these two modes can be explained 
entirely by examining their fundamental QNM frequen- 
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cies: CO 220 — ^210 = 0.08/Mf, a much slower frequency 
than CJ 330 - CJ 220 = 0.31/Mf. Prom Eqs. (59a), the final 
kick should be roughly ll-iccWch^ 77 / 4 / (^ 210 -^ 220 )! ~ 2 
times larger than at the end of inspiral, as seen in Fig. 12. 
It is somewhat more complicated to calculate the ampli- 
tude of the antikick for the other runs, which have mul- 
tiple dominant flux-pairs, and thus requires the summa- 
tion of Eqs. (59a)-(59c). For the NE+ 3 run we estimate 
^finai/^insp ~ 0.7 and for the NE§5 3 run t/finai/^insp ~ 0.85, 
in qualitative agreement with the numerical simulations. 

In general, we find the relative magnitude of the anti- 
kick is mostly dependent on the relative magnitudes of 
the S 21 and I 33 moments. When S 21 dominates, the 
ringdown rotation is slow and there is a small unkick, 
whereas a dominant 7 33 mode gives a rapidly rotating 
ringdown flux and thus a large unkick. Furthermore, 
from Eq. (51a), we see that for non-spinning BHs, both 
the S 21 and 7 33 modes share the same mass and fre- 
quency scaling, so the relative size of the unkick should 
be roughly independent of mass ratio. 

We want to get a more quantitative picture of how 
these flux vectors add constructively and destructively 
to give the total recoil velocity to support the analytic 
estimates presented above. Using v = f Fdt we can write 

f M = S (< " V) = *' F ’ (60) 

where v • v = 1. Breaking as above F up into the contri- 
butions of the dominant modes and then integrating in 
time gives 

^1,22 = yV F 21.22 dt 

v 22 ’ 33 = Jv-F 22 ’ 33 dt 

v 33,44 = y*. F 33,44 ^ (61) 

which add linearly to give to total recoil velocity: 

| v | = u 21 ’ 22 + v 22 ’ 33 4 - u 33>44 . (62) 

These different velocities are plotted in Figs. 13, with the 
same color scheme as in Fig. 12, along with the total in- 
tegrated recoil velocity in solid black curves. Also shown 
in Fig. 13 is the velocity t/ 32,33 (dashed blue curves), de- 
fined analogously to Eq. (61) for the S 32 / 33 * flux terms. 
The remarkably small contribution from this mode pair 
further justifies our focus on the more dominant pairs of 
Eq. (23) and Fig. 3. 

In the NE 2: + run, where the modal analysis shows 
the F 21}22 and F 33,44 flux terms canceling out, we see 


that the total integrated recoil velocity (black curves in 
Fig. 13) is almost entirely dominated by the F 22 ’ 33 flux 
(red curves). On the other hand, for the NE+ 3 run, the 
F 21 ’ 22 flux is much stronger, adding constructively with 
the F 22 ’ 33 flux. This has the effect of both increasing the 
peak velocity and also decreasing the relative strength 
of the unkick, due to the slow rotation frequency of the 
F 21,22 flux during ringdown, as described above. As ex- 
pected, the NE 2 5 3 run displays behavior intermediate be- 
tween these two extremes. The EQ_| run, however, is 

entirely dominated by the F 21,22 flux, and thus experi- 
ences almost no unkick, but rather drifts off slowly in 
a nearly constant direction, as seen in the bottom-right 
panel of Fig. 12. 


D. Application to non-planar kicks 

One of the most remarkable results from the recent 
renaissance in numerical relativity was the prediction of 
extremely large kicks from equal-mass BHs with spins 
pointing opposite to each other and normal to the or- 
bital angular momentum, producing a recoil out of the 
orbital plane [13-15]. While this configuration can pro- 
duce recoils of nearly 4000 km/sec, the analogous non- 
precessing configuation (EQ+_ run in this paper) give 
kicks of only ~ 500 km/sec in the case of maximal spin 
[10-12]. It turns out that the multipole analysis tools de- 
veloped above are ideal for understanding and explaining 
this remarkable difference. 

First, we should note that leading-order estimates of 
the linear momentum flux during inspiral suggest that 
the discrepancy should be less than a factor of two. For 
example, Eq. (3.31b) of Kidder [35] gives the spin-orbit 
contribution to the momentum flux for circular, Keple- 
rian orbits as 

Fso = ^ 2 ^3[n x A + (n x v)(v- A)], (63) 

with n and v being the normalized separation and veloc- 
ity vectors, respectively. For spins parallel to the angu- 
lar momentum, the term in brackets has magnitude A 2 , 
while for planar spins, it is 2|A| sin</>A, where <j > a is the 
angle between A and n. 

Not surprisingly, we see the exact same scaling from 
the multipole analysis of Eqs. (19), (20), (49), and one 
new moment: 

Sfo = -U\j^-r]Ru) 3 e~ i(l> (A x - iA y ). (64) 
Combining these equations, we get 
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FIG. 13: Relative contributions to the integrated recoil velocity from the different multipole mode-pairs. J 22 / 33 * ( re d curve) 
is the dominant mode for unequal-mass binary systems, while S 21 I 22 * (blue curve) dominates for spinning, equal-mass binary 
systems. Also plotted is the contribution from the S 32 / 33 * flux terms (blue dashed curve), demonstrating its very small 
contribution to the total recoil velocity. The data refer to the NE^i (upper left panel), NEi : ^_ (upper right panel), NEqo 3 
(lower left panel), and EQ+_ (lower right panel) runs. We denote with £ pea k the time at which I 22 reaches its maximum. 


and 


F z « g-E[— 28 S(/ 22 S 22 *)] = R 3 ( A * sin <fi - A v cos <j>). (66) 

f ~ 


So in both paradigms, we see that, when maximizing 
over sin 0 a, the planar-spin orientation should result in 
a recoil twice as large as the parallel-spin case, leaving a 
factor of roughly 4 difference unexplained. 

From Eqs. (65), (66) we see that the only relevant 
modes involved should be J 22 , S' 21 , and S 22 . In the left 
panel of Fig. 14 we plot the amplitude of 1 22 from the 
EQ^ — simulation, along with that of a planar-spin sim- 
ulation EQpianar- All other binary parameters and the 
initial conditions are the same. Remarkably, the mass 
quadrupole moments are nearly identical in both runs, 
suggesting that the overall dynamics of the inspiral and 
merger are the same. This is in fact quite reasonable 
since the total spin of the system is zero in both cases. 
However, we see in the right-hand panel of Fig. 14 that 
the peak amplitude of the S 22 mode is a factor of ~ 2.5 


greater than that of the S 21 mode from the EQ p i an ar and 
EQ + _ runs, respectively. 

Yet Eqs.’ (49), (64) suggest that these two modes should 
have the same magnitudes, at least during the inspi- 
ral phase, and presumably during the RD as well, since 
the RD amplitudes are completely determined by the 
mode amplitudes at the matching point. It appears from 
Fig. 14 that S 22 and S 21 do in fact have the same am- 
plitude at early times, but the relatively noisy data and 
short duration of the simulations make it impossible to 
say for certain. If this is the case, our best explana- 
tion for the sudden remarkable increase in the ampli- 
tude of S 22 is due to mode-mixing with the dominant 
1 22 mode, much like that of S 32 and 1 22 described above 
in Sec. IV B. This effect is apparently only important 
between modes with the same m- number, and may pos- 










FIG. 14: left panel: Comparison of the multipole amplitudes 1 22 for the two different equal- mass simulations: EQ p i an ar (solid 
line) and EQ+_ (dashed line), right panel: The S 22 amplitude from the planar-spins run (EQ p i anar , solid line) and the 
parallel-spins run (EQ+_, dashed line). 
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FIG. 15: left panel: Comparison of the linear momentum flux for the two different equal- mass simulations: EQ p i ana r (solid 
line) and EQ+_ (dashed line), right panel: The total recoil velocity from the planar-spins run (EQ p i anar , solid line) and the 
parallel-spins run (EQ+_, dashed line). 


sibly be explained by the fact that the QNMs are really 
expressed as spheroidal , not spherical harmonics [65, 66]. 

Lastly, from the ringdown contribution to the veloc- 
ity [Eqs. (58), (59)], we can understand another differ- 
ence between the planar- and parallel-spin orientations. 
Instead of having two different RD frequencies <7210 and 
0220 combine to give a slowly rotating flux vector, for the 
planar-spin case, we have two identical RD frequencies, 
giving precisely zero rotation to the RD flux. Further- 
more, as the spin vector A is precessing faster and faster 
in a positive direction around the orbital angular mo- 
mentum vector, even during the inspiral the two modes 
1 22 and S 22 become nearly locked in phase, producing 
a relatively long-duration burst of linear momentum flux 
in a single direction during the merger phase. Combined, 
these effects essentially straighten out the spiral curve in 
the lower-right panel of Fig. 12, providing another factor 


of ~ 1.5 of increased recoil velocity for planar spins. 

In Fig. 15 we show the combination of the above ef- 
fects. In the left panel, we plot the linear momentum 
flux from Eqs. (65), (66), showing the factor of two in- 
crease predicted by the Kidder formula, along with the 
factor of 2.5 increase in the amplitude of S 22 from mode- 
mixing. In the right panel, we plot the integrated recoil 
velocity for both runs, which includes the effect of flux 
rotation during the merger and inspiral phases, account- 
ing for another factor of ~ 1.6, giving a total discrepancy 
of u(EQ planar )/u(EQ + _) * 2.5 x 2 x 1.6 = 8. 

VII. DISCUSSION 

In this paper we analysed several numerical simula- 
tions of binary BH coalescence focusing on the physics of 
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the recoil. We developed tools, based on the multipolar 
expansion [41-45], that can be used as a diagnostic of the 
numerical results, and understand how the recoil velocity 
evolves during the inspiral, merger, and ringdown phases 
of the coalescence. 

We wrote explicit expressions for the linear momen- 
tum flux expressed in terms of radiative multipole mo- 
ments through l = 4, valid for generic spinning, pre- 
cessing BH binary systems. We show that these for- 
mulae are sufficient to obtain the final recoil with high 
accuracy. By comparing the amplitudes of the differ- 
ent multipole moments, we found that only three pairs 
of modes contribute significantly to most of the recoil, 
notably S 21 / 22 *, I 22 / 33 * and J 33 / 44 *. Those modes ac- 
count for the total recoil with an accuracy on the order 
of ~ 5 — 10% throughout the simulations, (see Figs. 8, 
9). 

The way in which the contribution from these three 
pairs of modes builds up is not trivial, since not only the 
relative amplitudes, but especially, the relative phases are 
also quite important. We find that the relative phases be- 
tween the three pair modes are nearly constant during the 
inspiral phase, but start evolving at the onset of the tran- 
sition from inspiral to RD (see Fig. 11). This late-time 
evolution can be described reasonably well with analytic 
formula obtained expressing the mode-pairs in terms of 
fundamental QNMs of a Kerr BH. We showed that it is 
the relative magnitude of the current-quadrupole mode 
S 21 and mass-octupole mode I 33 that determines the dif- 
ference between the recoil at the peak of the linear mo- 
mentum flux, and the final recoil velocity, i.e. , the amount 
of anti- kick. 


With the final goal of improving analytic PN mod- 
els, we also explored whether simple modifications of the 
Newtonian formula for the linear momentum flux, allow 
to match the numerical results all along the binary evo- 
lution. We found that, if we treat the binary radial sep- 
aration in the Newtonian multipole modes (41a)-(41e) 
with an effective radius, which is computed from the nu- 
merical simulations assuming that each multipole mode 
is described by a dominant frequency (see Fig. 1 ), the 
leading Newtonian modes reproduce quite well the nu- 
merical ones (see Figs. 6, 7) up to the end of the inspiral 
phase. We also found, confirming the results in Ref. [66], 
that a superposition of three QNMs, can fit very well the 
numerical waveforms beyond the peak of the radiation, 
all along the RD phase. 

The tools developed in this paper will be employed 
to improve the PN analytic models [27] and its siblings, 
notably the EOB model [28, 29, 31], so that they can 
accurately interpolate between numerical results. In this 
way, fast Monte Carlo simulations will be able to predict 
recoil distributions from BH mergers with uncertainties 
smaller than in Ref. [23]. Those recoil distribution can 
be included in simulations of hierarchical merger models 
of supermassive BHs providing more robust predictions. 
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